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Considerable progress has been recently made in the theoretical understanding of the colossal 
magnetoresistance (CMR) effect in manganites. The existence of inhomogeneous states has been 
shown to be directly related with this phenomenon, both in theoretical studies and experimen- 
tal investigations. The analysis of simple models with two competing states and a resistor network 
approximation to calculate conductances has confirmed that CMR effects can be theoretically repro- 
duced using non-uniform clustered states. However, a direct computational study of the CMR effect 
in realistic models has been difficult since large clusters are needed for this purpose. In this paper, 
the recently proposed Truncated Polynomial Expansion method (TPEM) for spin-fermion systems 
is tested using the double-exchange one-band, with finite Hund coupling Jh, and two-band, with 
infinite Jh, models. Two dimensional lattices as large as 48x48 are studied, far larger than those 
that can be handled with standard exact diagonalization (D1AG) techniques for the fermionic sec- 
tor. The clean limit (i.e. without quenched disorder) is here analyzed in detail. Phase diagrams are 
obtained, showing first-order transitions separating ferromagnetic metallic from insulating states. A 
huge magnetoresistance is found at low temperatures by including small magnetic fields, in excellent 
agreement with experiments. However, at temperatures above the Curie transition the effect is much 
smaller confirming that the standard finite-temperature CMR phenomenon cannot be understood 
using homogeneous states. By comparing results between the two methods, TPEM and DIAG, on 
small lattices, and by analyzing the systematic behavior with increasing cluster sizes, it is concluded 
that the TPEM is accurate to handle realistic manganite models on large systems. Our results 
pave the way to a frontal computational attack of the colossal magnetoresistance phenomenon using 
double-exchange like models, on large clusters, and including quenched disorder. 



PACS numbers: 75.47.Lx, 75.30.Mb, 75.30.Kz 

I. INTRODUCTION 

The study of transition metal oxides (TMOs) has 
been among the most important areas of investigations 
in condensed matter physics in the last two decades. 
The excitement in this context started with the high- 
temperature superconductors and was later followed 
by the discovery of the colossal magnetoresistance in 
manganites as well as a plethora of equally interesting 
phenomena in several other oxides. Strong correlations 
(i.e. large electron-electron or electron-phonon couplings, 
or both) play a major role in the physics of these materi- 
als. In addition, the presence of nearly degenerate states 
renders some of these oxides highly susceptible to ex- 
ternal perturbations. In fact, TMOs appear to share a 
phenomenology similar to that of complex systems, with 
nonlinearities and sensitivity to details^ 

We focus in this work on the manganites, area where in 
recent years considerable progress has been made, both 
in theory and experiments^ In the late 90's, it was pre- 
dicted that many Mn-oxides should be inhomogeneous 
at the nanoscale, due to the unveiling of tendencies to- 
ward electronic phase separation^ On the experimental 
front, the evidence for inhomogeneous states was rapidly 



gathered and it is by now widely accepted, with build- 
ing blocks typically having small nanoscale sizes£ Sub- 
sequent theoretical work showed that similar tendencies 
can also occur after the inclusion of quenched disor- 
der effects - caused for instance by chemical doping- 
near first-order phase transitions^ The key influence of 
quenched disorder was also observed in simulations of 
the one-band model for manganites including coopera- 
tive phonons^ and also for two bands with Jahn- Teller 
phonons^> This key role of quenched disorder postulated 
by theory was observed experimentally using a Mn-oxide 
compound that can be prepared both in crystal ordered 
and disordered forms^S*ii Remarkably, only the latter 
presents a CMR effect. 

While the presence of quenched disorder was theoreti- 
cally found to generate metal insulator transitions simi- 
lar to those found experimentally, the actual existence of 
large magnetotransport effects is difficult to test in un- 
biased theoretical studies. Using toy models that have 
phase competition and first-order transitions, and sup- 
plementing the investigations with a random-resistor net- 
work approximation, huge magnetoresistances were ob- 
tained in resistance vs. temperature profiles in excellent 
agreement with experiments. 12 However, it is certainly 
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desirable to carry out similar investigations in more re- 
alistic models, of the double-exchange variety, and with 
quantum mechanical procedures to calculate the conduc- 
tances. Alas, this task is tremendously difficult with 
standard computational methods that rely on the exact 
diagonalization of the fcrmionic sector and the Monte 
Carlo simulation of the classical i 2 g spinsi^ The effort 
in this context grows like iV 4 , with N the number of 
sites, severely limiting the clusters that can be stud- 
ied. Since theory suggests that physics related with per- 
colative phenomena is anticipated upon the introduction 
of magnetic fields in nano-clustered states (2*^ large clus- 
ters must be used for accurate simulations. Fortunately, 
important progress has been made in recent years to- 
ward the development of a novel method to carry out 
these investigations. 13 ! 14 ' 15 ! 16 The technique, TPEM, has 
a CPU time that scales like N and, as a consequence, it 
is a promising technique for these investigations. Previ- 
ous studies have shown that the Jjj=oo one-band model 
with and without quenched disorder can be accurately 
studiedii 3 ^ 1 ^ 5 ^* 1 ^ However, the method has not been 
tested yet under some of the more severe circumstances 
needed for a realistic study, namely the use of two active 
orbitals per site (i.e. the doubly degenerate e g sector of 
Mn ions) and/or with a finite Hund coupling. 

It is important to remark that there are at least two 
other independent methods that have been proposed to 
improve on the exact diagonalization of the fcrmionic sec- 
tor approach: (1) The hybrid Monte Carlo technique of 
Alonso et a/., 18 inspired in lattice gauge theory, that com- 
fortably allows calculations on lattices up to 10 3 sites in 
the limit of Jh=oo and for one orbital^ and with a linear 
with TV increase in effort, and (2) the method of Kumar 
and Majumdar— that has been already applied to a va- 
riety of models reaching 32x32 clusters for one orbital 
and at Jh=oo. The scaling of this method is TV 3 . Our 
choice of the TPEM is motivated by the perception that 
a linear with N cost is needed to handle the anticipated 
percolative physics that emerges when quenched disor- 
der and phase competition occurs. Also it seems easier 
to implement than method (1) where auxiliary fields are 
needed. Nevertheless, this is not a critique on methods 
(I) and (2): they should be strongly pursued as well. 
Only future work can decide which method is the best 
for the type of problems described here. 

It is the main purpose of this paper to present a 
systematic study of the TPEM applied to models that 
are widely believed to be realistic for manganite inves- 
tigations, in the clean limit. The conclusions indicate 
that the technique works properly for these investiga- 
tions, opening an avenue of research toward the ultimate 
goal of conducting a fully realistic simulation of the two- 
band double-exchange model including quenched disor- 
der. The present results include a detailed analysis of 
both metallic and insulating phases on lattices as large as 
48x48 sites, about 20 times larger in number of sites than 
it is possible to handle with exact diagonalization tech- 
niques. In addition, here it is discussed the existence of 



a huge magnetoresistance at low temperatures, in agree- 
ment with experiments for (Ndi_j / Smy) 1 / 2 Sr 1 / 2 Mn03. 21 
This phenomenon was theoretically studied before by 
Aliaga et al~ as well, although on much smaller systems. 
Nevertheless, the conclusions are similar: this somewhat 
exotic "low temperature" CMR phenomenon can be un- 
derstood as a natural consequence of the existence of 
a first-order metal-insulator transition in the phase di- 
agram and, thus, a clean-limit simulation is sufficient for 
this purpose. However, in our investigations it is also con- 
firmed that these clean limit simulations are not able to 
generate the standard CMR effect above the Curie tem- 
perature using states that are homogeneous. Quenched 
disorder or strain fields are likely important for this pur- 
pose, which will be the subject of near future efforts. 

The organization of the paper is simple. The theo- 
retical aspects of the TPEM are briefly reviewed in Sec- 
tion [H] The systematics of the TPEM in the case of the 
one-band model, with emphasis on the dependence with 
parameters and size effects is discussed in Section 11111 
We also present physical results related with large mag- 
netoresistance effects found at low temperatures, in large 
clusters. Then, in Section Hvl the emphasis shifts to the 
two bands model, with an analogous study of TPEM per- 
formance and effects of magnetic fields. Conductances 
are evaluated for both models using the TPEM and rea- 
sonable results are observed. Overall, it is concluded that 
the TPEM is adequate for a frontal future attack of the 
CMR phenomenon using realistic models and quenched 
disorder. 



II. REVIEW OF THE TPEM 

For completeness, here a brief review of the TPEM 
is presented following closely Ref. Consider 
a model defined by a general Hamiltonian, H = 
Ylijap c \a H ia,ji3(4>)c-]P, where the indices i and j denote 
a spatial index, while a and /3 are internal degree(s) of 
freedom, e. g. spin and/or orbital. As in the case of 
spin-fermion models, the Hamiltonian matrix depends on 
the configuration of one or more classical fields, repre- 
sented by (j). Although no explicit indices will be used, 
the field(s) are supposed to have a spatial dependence. 
The partition function for this Hamiltonian is schemati- 
cally given by: Z = f d^J2 n ( n \ ex P(-/ 3i ^('/') + P^N)\n) 
where \n) are the eigenvectors of the one-electron sector, 
and the <f> integral denotes the integration over all the 
classical fields. Here 0=l/{kB,T) is the inverse tempera- 
ture. The number of particles (operator N) is adjusted 
via the chemical potential fi. 

The procedure to calculate observables (energy, den- 
sity, action, etc) is the following. First, an arbitrary 
configuration of classical fields <f> is selected as a starting 
point. The Boltzmann weight or action of that configura- 
tion, S ((/>), is calculated by diagonalizing the one-electron 
sector. Then, a small local change of the field configura- 
tion is proposed, so that the new configuration is denoted 
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by <j>' and its action is recalculated to obtain the differ- 
ence in action AS — S(<p') — S(<f>). Finally, this new 
configuration is accepted or rejected based on a Monte 
Carlo algorithm, such as Metropolis or heat bath, and 
the cycle starts again. In summary, in the standard algo- 
rithm (that we will denote here as DIAG) the observables 
are calculated using an exact diagonalization of the one- 
electron sector for every classical field configuration, and 
Monte Carlo integration for those classical fields^ 

The TPEM replaces the diagonalization for a polyno- 
mial expansion as discussed below (details can be found 
m Refs-EHHl. It will be assumed that the Hamilto- 
nian H (0) is "normalized" , which simply implies a rescal- 
ing in such a way that the normalized Hamiltonian has 
eigenvalues e v £ [—1,1]. Simple observables can be di- 
vided in two categories: (i) those that do not depend 
directly on fermionic operators, e. g. the magnetization, 
susceptibility and classical spin-spin correlations in the 
double exchange model, and (ii) those for which a func- 
tion F(x) exists such that they can be written as: A(4>) — 
f , F(x)D(4>, x)dx, where D(<p, e) = Y^ v S{t{4>) — e v)i and 
e„ are the eigenvalues of H{4>), i. e. D(4>, x) is the density 
of states of the system. 

For category (i) , the calculation is straightforward and 
simply involves the average over Monte Carlo configura- 
tions. Category (ii) includes the effective action or gener- 
alized Boltzmann weight and this quantity is particularly 
important because it is calculated at every MC step to 
integrate the classical fields. Therefore, we first briefly 
review how to deal with this type of observables. As 
Furukawa et al. showed, it is convenient to expand the 
function F(x) in terms of Chebyshev polynomials in the 
following way: F(x) = J2m=o fm,T m (x), where T m (x) is 
the m— th Chebyshev polynomial evaluated at x. Let 
a m = 2 — (5-m.o- The coefficients f m are calculated with 
the formula: f m = a m F(x)T m (x) / (n\/l — x 2 ). The 
moments of the density of states are defined by: 

N dim 

fi m (<t>) = J2 {v\T m (H(4>))\v), (1) 

where Ndim is the dimension of the one-electron sector. 
Then, the observable corresponding to the function F(x), 
can be calculated using 

A{cj>) = Y d fm^m{4>)- (2) 
m 

In practice, the sum in Eq. J2Jl is performed up to a 
certain cutoff value M, without appreciable loss in ac- 
curacy as described in Refs. ^| and Q, and as exten- 
sively tested in the main sections of this paper for re- 
alistic manganite Hamiltonians. The calculation of /i m 
is carried out recursively using \v;m) = T m {H)\v) = 
2H\v; in — 1) — \v; m — 2). These same vectors are used 
to calculate the moments. The process involves a sparse 
matrix- vector product, e. g. in T m (H)\v) , yielding a cost 
of 0(N 2 ) for each configuration, i. e. 0(N 3 ) for each 



Monte Carlo step. This represents an improvement in a 
factor N compared with DIAG. 

In addition, an extra improvement of the method de- 
scribed thus far has been proposed^ based on two trun- 
cations: (i) of the sparse matrix-vector product and (ii) 
of the difference in action for local Monte Carlo updates. 
The resulting algorithm has an expected CPU time grow- 
ing like N. The first truncation is possible because of 
the form of the vectors \u; m). Indeed, for m=0 these are 
simply basis vectors that can be chosen with only one 
non-zero component. The m=l vector is obtained by ap- 
plying Ti(H)=H to the basis vector. Since H is sparse 
(consider for instance a nearest-neighbor hopping), the 
vector \v; m = 1) will have non-zero elements only in the 
vicinity of v. For general m, the non-zero elements will 
propagate in what resembles a diffusion process. Note 
that we only have to keep track of the non-zero indices 
to perform the sparse matrix- vector product. Since we 
are only discarding null terms, this truncation does not 
introduce any approximation. It is possible to go a step 
further and discard not only zero elements but elements 
smaller than a certain threshold denoted by e pr . In this 
case, the results are approximate, but the exact results 
are recovered in the limit e pr — > 0. In this paper, the 
dependence of results with various values for this cutoff 
(and the one described below) is discussed. 

The second truncation involves the difference in effec- 
tive action, which is calculated very frequently in the 
Monte Carlo integration procedure. The function cor- 
responding to the effective action for a configuration 4> 
is defined by F s (x) = — log(l + exp(— j3{x — fi))) and 
S((f>) admits an expansion as in Eq. J5J, with coefficients 
corresponding to F s (x). In practice, only the dif- 
ference in action, AS=S(<fi') — S(4>) has to be computed 
for every change of classical fields. Since this operation 
requires calculating two sets of moments, for (f> and </>', 
the authors of Ref . ^| have also developed a truncation 
procedure for this trace operation controlled by a small 
parameter, et r - This truncation is based on the observa- 
tion that if cf> and </>' differ only in very few sites (as is 
the case with local Monte Carlo updates), then the cor- 
responding moments will differ only for certain indices 
allowing for a significant reduction of the computational 
effort. Again, the exact results are recovered for e tr — > 0, 
so this approximation is controllable. 

Another key advantage of the TPEM is that it can be 
parallelized, because the calculation of the moments in 
Eq. Q for each basis ket \v) is independent. Thus, the 
basis can be partitioned in such a way that each proces- 
sor calculates the moments corresponding to a portion 
of the basis. It is important to remark that the data to 
be moved between different nodes are small compared to 
calculations in each node: communication among nodes 
is mainly done here to add up all the moments. The pos- 
sibility of parallelizing this algorithm can be contrasted 
with the conventional method where a matrix diagonal- 
ization is performed at every Monte Carlo step; in that 
case the calculations must be serial because it is difficult 
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to make an efficient parallelization of the matrix diago- 
nalization. 



III. RESULTS FOR THE ONE BAND MODEL 
A. Definition 

After introducing the computational method, we will 
now focus on its performance starting with the one- 
band model for manganites. Historically, this model was 
among the first proposed for Mn-oxides and it is still 
widely used, although it does not incorporate the two 
orbitals e g of relevance in Mn ions. This more involved 
two-band version will be studied in the next section. 

The one band Hamiltonian used in this study is given 
by: 

Hu = -t (4,a c j> + h.c.) -JnYl c\ a a a ^a^ 

(ij),a i,a,P 

+ Jaf ^ S l • Sj, (3) 

m 

where cj creates an electron at site i with spin a, a 
are the Pauli spin matrices, Si is the total spin of the 
t 2g electrons (assumed to be localized and classical), (ij) 
indicates summing over nearest neighbor sites, t is the 
nearest neighbor hopping amplitude for the movement of 
electrons (t also sets the energy unit), Jh > is the Hund 
coupling, and Jaf > is the antiferromagnetic coupling 
between the localized spins. The study carried out in 
this manuscript is based on the clean limit, namely the 
couplings that appear in the Hamiltonian do not have a 
site index, which would be necessary if quenched disorder 
is incorporated. The study of realistic models including 
disorder will be carried out in a future investigation, since 
it represents an order of magnitude extra effort due to the 
average over disorder configurations. 

Here and in the rest of the paper, spatial indices will 
always be denoted without arrows or bold letters inde- 
pendently of the dimension. Also the notation i + j is 
meant to represent the lattice site given by the vectorial 
sum of the vectors corresponding to i and j, respectively. 

B. TPEM performance 

1. Test of the TPEM in Small Systems 

As explained before, the TPEM has three controlling 
parameters: M, e pr , and et r - In the limit when the 
first parameter runs to infinity, and the other two to 
zero, the exact results are recovered. For the TPEM to 
be useful, accurate results must be obtained for values 
for these parameters that allow for a realistic computa- 
tional study. In Fig. ^a), the dependence of the zero- 
momentum spin structure factor, of relevance for ferro- 
magnetism, is shown vs. temperature, using a 12 2 cluster 
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FIG. 1: (color online). Dependence of the TPEM algorithm 
results on the number of terms M in the expansion. Shown 
are the spin structure factors at the momenta characteristic 
of (a) the FM state and (b) the Flux state, normalized to 
1 and 0.5, respectively, for the perfect states. Results are 
obtained on a 12 x 12 lattice and they are compared with the 
numbers gathered using the exact diagonalization algorithm, 
all calculated at density (n) — 0.5. (a),(c) correspond to 
Jaf ~ 0.0 and (b),(d) to Jaf = 0.1. Measurements were 
taken every 10 steps of a MC run of 2000 total iterations, after 
discarding 2000 steps for thermalization. A random starting 
configuration is used for each T. In (a), at the most difficult 
temperature, T=0.06, where critical fluctuations are strong, 
the result shown was confirmed using several different starting 
configurations, including ordered ones. The average S(0, 0) 
obtained by this procedure was very similar among the several 
starts. The shown error bars at this temperature and M=30 
mainly arise from the expected critical fluctuations. In (b), 
a good convergence at T=0.04 is only achieved by using 40 
moments with e pr = 10 -7 and etr = 10 -5 and the result is 
shown with an orange star just below the exact result. In 
(c), the TPEM parameters used are M = 30, e pr = 10~ 5 
and etr = 10 -3 , except at T=0.06, where the convergence is 
achieved by using 40 moments with e pr = 10 -7 and etr = 
10 _J . In (d), all the results shown were obtained with M = 
30, e pr = 10~ 5 and etr = 10~ 3 , and solid lines represent the 
DIAG results while the dashed lines are the TPEM results. 



and the values of Jh and Jaf indicated. In this case, it 
is expected that a FM state will form at low tempera- 
tures, as observed numerically Results for many values 
of M are shown, at fixed values of e pr , and etr- Clearly, 
M=10 only captures the low and high temperature lim- 
its, but it is not accurate near the critical temperature. 
The results for M=20 are much better, but still there is a 
visible discrepancy near the region where S(0, 0) changes 
the fastest. However, for M=30 and 40, fairly accurate 
results are obtained. In Fig. [He), it is shown that even 
the spin correlations at the largest distances are accu- 
rately reproduced with M =40 terms in the expansion. 

In Fig. [IJb), similar results are presented but now 
for the spin-structure-factor corresponding to the "Flux" 
phase - nearest-neighbor spins at 90 degrees forming a 
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FIG. 2: (color online), e dependence of the TPEM algorithm 
results (M=30) for a 12 x 12 lattice, compared with the exact 
(DIAG) results. Both are obtained at density (n)=0.5, using 
(a) Jaf = 0.0 and (b) Jaf = 0.1. Measurements were taken 
at every 10 steps of a 2000 MC steps total run, after 2000 
steps for thermalization. In (b), the convergence at T=0.04 
is achieved by using 40 moments with e pr = 10" 7 and etr = 
10 -5 , as shown with an orange star just below the exact result. 
The starting configurations and error bar convention is as in 
FigE] 



staggered arrangement of nonzero plaquette fluxes. This 
phase appears at n=0.5 with increasing Jaf, as reported 
in previous investigations, 22 In this case, M=10 and even 
20 produces results dramatically different from those of 
DIAG. However, for M=30 discrepancies are observed 
only in a range of temperatures near the critical tran- 
sition, with low and high temperatures under control. 
Finally, M=40 leads to very accurate results, as in case 
(a). Even the spin correlations are under well control for 
this number of terms in the expansion (see Fig. ^d)). 
The range M~30-40 appears systematically in our in- 
vestigations, and it is expected to provide safe values of 
M for studies of the type of spin-fermion models under 
investigation in manganites. 

In Fig. |3 a study of the dependence of results with 
e pr and etr is presented, working at M=30. From Fig. [21 
clearly there are large e's that lead to unphysical results, 
but with decreasing values an accurate evaluation of ob- 
servables is reached. In this and other investigations, val- 
ues such as e pr =10~ 5 and e tr =10~ 3 are generally found 
to be accurate, with only a few exceptions. 



2. Dependence of Results on Lattice Sizes 

An approximate method that depends on some param- 
eters, such as in the case of the TPEM, is practical only 
if by fixing those parameters on small systems, their val- 
ues still provide accurate numbers as the lattice sizes in- 
crease. A qualitative way to carry out this test is to 
perform the studies on large clusters and see that all the 
trends and approximate numbers remain close to those 




FIG. 3: (color online). Lattice size dependence of the TPEM 
algorithm results for the lattices and parameters shown, work- 
ing at density (n)=0.5 and using (a) Jaf = 0.0 and (b) 
Jaf = 0.1. Measurements were taken at every 10 steps of a 
2000 MC total steps run, after discarding 2000 MC steps for 
thermalization. For the 20x20 lattice and larger, the start- 
ing configuration used was a perfect FM state. In (b), the 
starting configuration used is a random one except at T=0.04 
where the starting configuration is chosen to be a perfect flux 
state for faster convergence. In addition, for lattices except 
40 x 40, the convergence at T=0.04 is achieved by using 40 
moments with e pr = 10~ 7 and etr = 10~ 5 , while for other 
temperatures M = 30 with e pr = 10~ 5 and etr = 10~ 3 were 
sufficient. For the lattice 40 x 40, M = 40 with e pr = 10~ 6 
and etr = 10~ 4 were used for all temperatures. 



known to be accurately obtained on small systems, or 
expected from other techniques or physical argumenta- 
tions. Figure Efa,b) supports the notion that TPEM in- 
deed behaves properly in this respect, namely the range 
of M and e's identified in the previous subsection are suf- 
ficient to produce qualitatively similar results even when 
the number of sites grows by a factor 10. In (a), the ex- 
pected size dependence corresponding to a second order 
FM transition is found. For a 40x40 cluster, the Curie 
temperature appears located at T~0.07. In (b), the size 
dependence is almost negligible. The transition is far 
sharper for the paramagnetic-flux transition, as already 
noticed in Fig. EJb) . This is an intriguing feature that 
will be investigated in future work: while the first-order 
low-temperature metal-insulator transitions are clear and 
well established in realistic models for manganites, the 
presence of first-order transitions between ordered and 
disordered phases varying temperature is far less obvious, 
and TPEM studies on large lattices can properly address 
this issue. For our current purposes, here it is sufficient 
to state that the TPEM appears to behave properly with 
increasing lattice sizes, both in metallic and insulating 
regimes. 

At this point a clarification is important. In principle, 
two-dimensional systems should not show true critical 
temperatures due to the Mermin- Wagner theorem. How- 
ever, it is well known that in systems where the Mermin- 
Wagner theorem applies, such as the two-dimensional NN 



TABLE I: Comparison of the CPU times for the algorithms 
indicated, using an Intel Pentium 4 (clock speed 3.06Gz) com- 
puter. Shown are results for different square lattices of size 
L x L, assuming 2000 MC steps for thermalization, and 2000 
MC steps for measurements (taken every 10 MC steps) . Since 
the TPEM can be parallelized, some results were obtained us- 
ing more than one CPU, as indicated. 
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Heisenberg model, the antiferromagnetic spin correla- 
tions exponentially diverge with decreasing temperature. 
An exponential behavior, defines via the exponent a tem- 
perature scale T* below which the correlations are much 
larger than any lattice size that can be practically stud- 
ied numerically. This may seem like a problem, but it is 
not: very large correlation lengths also render the system 
very susceptible to small perturbations. In particular, we 
have shown that tiny deviations from the fully symmetric 
Heisenberg model, such as introducing Ising anisotropies, 
stabilize T* into a true critical temperature. In fact, 
simulations performed with Ising anisotropies typically 
reveal no important differences with the results obtained 
with fully vector models on finite systems. Small cou- 
plings in the third direction play a similar role. As a 
consequence, for all practical purposes the critical be- 
havior observed in the present studies describes properly 
the expected physics of manganite models, which are al- 
ways embedded in three dimensional environments, and 
that have small anisotropies. A final note on this sub- 
ject: The CE phase of manganites can show a finite- 
temperature transition even in two dimensions, since the 
order parameter for charge order can be Ising type. 

The CPU time needed to obtain the results shown in 
this subsection follows the expected trends reported in 
previous investigations (see Table [IJ. In particular, the 
TPEM time needed for a 32x32 cluster is comparable to 
the DIAG time on a 12x12 cluster, a very encouraging 
result. Of course, this comparison will be even more fa- 
vorable to the TPEM with increasing number of CPU's 
for parallclization. 



C. Phase Diagram 

Using the TPEM, the phase diagram of the one- 
band model for manganites at (n}—0.5 was obtained (see 
Fig. The transition between the FM and Flux states 
at low temperatures is of first order. In fact, in the ab- 
sence of quenched disorder the zero temperature result 
can be obtained by using the perfect classical spin config- 
urations for both the FM and Flux states, and calculating 
their energy vs Jaf (not shown). By this procedure the 
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FIG. 4: (color online). Spin structure factor vs. T for 
Jaf=0.05 using the lattice sizes and parameters shown in the 
figure. 
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FIG. 5: (color online). Phase diagram at (n)=0.5, varying 
temperature and Jaf- Results are shown for a 12 x 12 lattice 
using both DIAG and TPEM techniques, and for larger lat- 
tices using TPEM, as indicated. The origin of the tilting of 
the first-order low-temperature FM-Flux line is explained in 
the text. 



zero-temperature critical Jaf was found to be close to 
0.03. Raising the temperature, this transition line is not 
vertical, but has a tilting. Figure shows that the esti- 
mated critical temperatures do not present severe size ef- 
fects, and the TPEM can be comfortably used at least up 
to 40x40 clusters. The presence of a first-order transition 
in the competition between the FM and Flux phases is 
in qualitative agreement with several previous investiga- 
tions that have shown similar trends both for the one and 
two bands models, at any electronic density^ This tran- 
sition is expected to be severely affected by the influence 
of quenched disorder, and this issue will be investigated 
in the near future. 
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FIG. 6: (color online). Examples of the criteria used in the 
calculation of the critical temperatures in Fig. [£] (a) and 
(c): Spin structure factors at the momenta of relevance vs. 
temperatures, for many values of Jap, as indicated, using the 
DIAG technique on a 12 x 12 cluster, (b) and (d): Same as 
in (a) and (c) except the technique used here is the TPEM. 
Shown are some of the results for 12 x 12 and 20 x 20, as 
indicated in the figure. 



The critical temperatures in Fig.[3]were estimated from 
the behavior of the spin structure factors at the two mo- 
menta of relevance for the FM and Flux phases, as shown 
in Fig. El 
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FIG. 7: (color online), (a) Density-of-states calculated for a 
perfect Flux state using both the DIAG and TPEM methods 
for a lattice of size 12 x 12. In order to get the DOS accu- 
rately, even removing the Gibbs oscillations, one needs larger 
number of moments than what is usually required for other 
observables. (b) Monte Carlo results for the density-of-states 
obtained from simulations performed on a 40x40 lattice. In 
this case the last configuration of the MC run has been used 
to calculate the density-of-states at T=0.02. 



D. Density of States 

In this section, it is shown that the density-of-states 
(DOS) can be reproduced properly by the TPEM. This 
is nontrivial, since it may be suspected that a method 
based on an expansion of the DOS may have problems 
in an insulating phase due to the rapid changes in the 
DOS near the gap. To our knowledge, this is the first 
time that the TPEM is applied to an insulator. The 
discussion in this subsection shows that the technique 
works satisfactorily. 

In general, the density-of-states for a configuration of 
classical fields (j> is given by 

AVM = 5>(^'-6 A ), (4) 

A 

with uj' = {a} — b)/a, and where a and b are the pa- 
rameters that normalize the Hamiltonian in such a way 
that the new eigenvalues are in the interval [—1,1]. 
These constants are given by a=(E max — E min )/2 and 
b=(E max + E m i n )/2, where E max and E min are the max- 
imum and minimum eigenvalues of the Hamiltonian. 
Then, following the discussion of Section [DJ the corre- 
sponding function F(x) for the density-of-states in the 



expression 

A{4>)= J F(x)D(<P,x)dx (5) 
is the (^-function 

F(x) = 8(w'-x). (6) 

In the expansion f(i)=Em=otfm(4 b Y 
using the expression for the coefficients 

f m = f , a m F(x)T(x)/(TT\ / l — x 2 ), the final result 
for the density-of-states becomes 

N(u) = E m °mWMfl (7) 

7TV 1 — UJ 2 

This sum is truncated to a certain cutoff M and such an 
abrupt truncation results in unwanted Gibbs oscillations, 
as shown, e.g., in Fig.Qfor M=30. This problem may be 
avoided by multiplying the moments by dumping factors 
if needed^ 

The results for the density-of-states of the Flux phase 
are shown in Fig. [7| Clearly, even at M=30 there is a 
very good agreement between the DOS calculated exactly 
and with TPEM (with the exception of the in-gap Gibbs 
oscillations). Increasing M further, even this effect dis- 
appears. Using a 40x40 lattice, the results are almost 




FIG. 8: (a) Conductance and (b) resistance (1/conductance) 
vs. temperature for a 12 x 12 lattice, calculated with both 
DIAG and TPEM algorithms showing that the results agree. 
The couplings used are Jaf=0.0 and Jaf=0.1 as indicated, 
and the density is (n}=0.5. The convergence at T=0.04 is 
achieved by using 40 moments with £ pr =10 -7 and etr=10 -5 , 
as discussed in previous figures captions. 



FIG. 9: (color online), (a) Conductance and (b) resistance 
(1/conductance) vs. temperature calculated with the TPEM 
algorithm at Jap = 0.0 and Jaf = 0.1, (n)=0.5, and for 
the cluster sizes shown in the figure. The figure shows that 
the conductance does not suffer from strong size effects. The 
convergence at T=0.04 for the 12 x 12 lattice was achieved 
by using 40 moments with e pr = 10 -7 and etr = 10 -5 , as 
discussed elsewhere. 



the same as those observed on the smaller system. It is 
concluded that the TPEM can produce the DOS of insu- 
lating states accurately, and the method can be used to 
study phase competition between metals and insulators. 

E. Conductances: Comparison TPEM vs. DIAG, 
and Results with Increasing Lattice Sizes 

To compare theory with experiments, it is crucial to 
evaluate the conductance of the cluster under study. Its 
temperature and magnetic field dependence will clarify 
whether the double-exchange models for manganites con- 
tain the essence of the CMR phenomenon. The conduc- 
tance calculation here follows the steps previously exten- 
sively discussed by Verges et ai, and it basically relies 
on the Landauer formalism that links conductance with 
transmission. We refer the readers to original references 
for more details (see for instance Ref. |24}) . In Fig. [H] the 
conductance and its inverse (resistance) are shown as a 
function of temperature for the model on a 12x12 clus- 
ter that can be solved both with DIAG and TPEM. The 
agreement between the results obtained with both tech- 
niques is excellent, at the two values of Jaf shown (one 
in the FM and the other in the Flux phase). Thus, the 
conductance calculation does not present an obstacle in 
the use of the TPEM . 

With increasing lattice size, the TPEM conductance 
behaves smoothly and the finite-size effects are small 
(see Fig. |SJ, with the only exception of the insulating 
Flux phase regime at low temperatures where the 12x12 
cluster results appear appreciably different from those on 
larger systems. Considering the small value of the con- 
ductance in this insulating regime and the subsequent 
convergence of the Flux-phase resistance between the 
20x20 and 32x32 clusters, this appears to be only a mi- 



nor issue. 



F. Influence of Magnetic Fields in the Clean Limit 
and Partial Conclusions 

As discussed in the introduction, it is important to in- 
vestigate if the models studied here, in the clean limit, 
present a large magnetoresistance effect. Previous stud- 
ies by Aliaga et al& on 4x4 clusters, suggested that the 
"low temperature" large magnetoresistance experimen- 
tally observed in some manganites^ can be explained by 
double-exchange models in the clean limit. This result 
is important and deserves to be confirmed using larger 
clusters. Here, the case of the one-band model is ana- 
lyzed, with results shown in Fig. 1101 (two bands will be 
studied later in this paper). The value of Jaf was cho- 
sen to be on the insulating side (Flux phase) of the phase 
diagram Fig. El but close to the first-order transition sep- 
arating the metal from the insulator. The application of 
'small' magnetic fields favors the FM state over the Flux 
state and that manifests as a sharp transition from the 
metal to the insulator, for values of the magnetic field 
that appear abnormally small in the natural units of the 
problem. Thus, this model presents a huge negative mag- 
netoresistance, an encouraging result that shows theory 
is in the right track to understand manganites. The ef- 
fect shown in Fig. llOl is caused by the proximity in energy 
of two states with quite different properties, i.e. there is 
a hidden small energy scale in the problem. 

However, note that the standard large finite- 
temperature magnetoresistance traditionally studied in 
Mn-oxides cannot be understood with clean limit mod- 
els, as shown in Fig. the zero magnetic-field resis- 
tivity does not have the large peak near Curie temper- 
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FIG. 10: (color online), (a) Conductance vs. temperature 
for lattices of sizes 20 x 20 and 32 x 32, with and without 
magnetic fields, (b) Resistance vs. temperature calculated by 
taking the inverse of the conductance in (a) . Close to the first- 
order transition in the phase diagram, a small magnetic field 
can destabilize the insulating Flux state into a metallic FM 
state. For each temperature, 1000 thermalization and 2000 
measurements MC steps were used, with actual measurements 
taken at every 10 steps. 



atures characteristic of CMR manganites. Future work 
using the TPEM will analyze whether this more tradi- 
tional CMR effect can be obtained including quenched 
disorder. 

Overall, it can be safely concluded that the study of 
the one-band model with the TPEM has proven that the 
technique works properly, and that the FM vs. Flux 
competition occurs via a first-order metal-insulator tran- 
sition. This regime is ideal for the analysis of the influ- 
ence of quenched disorder in future calculations. 



IV. RESULTS FOR THE TWO-BAND MODEL 

A. Definition 

In this effort, the two-band model for manganites was 
also investigated using the TPEM. The two bands arise 
from the e g bands that are active at the Mn ions in Mn- 
oxides, as extensively discussed before^ The overall con- 
clusion of this section is that the TPEM is also a good 
approximation to carry out computational studies, con- 



clusion similar to that reached for only one active band. 
The Hamiltonian for this model is£ 



+ ^y^XQliPi + QljTxi + QztT zi ) 
i 

a=3 



(8) 



where the factor that renormalizes the hopping in the 
Jh=oo limit is 



S(9i,(tH,0j,^j) = cos(-^)cos(^-) 

+ sin(|)sin(|)e-^-^). (9) 



The parameters i" / are the hopping amplitudes between 
the orbitals 7 and 7' in the direction a. In this pa- 
per, we restrict ourselves to two dimensions, such that 
t x aa = -V3t* ab = -V3q a = 3t£, - 1, and tj, = ^ = 
V3t ba — 3t bb = 1. Qu, Qn and Q$i are normal modes of 
vibration that can be expressed in terms of the oxygen 
coordinate 



Qli = -^=[(u ii2 - Ui- ZtZ ) + (u^ x - Ui- X>x ) 
+ ( u i,y ~ u i-y,y)]> 



Q 2i 



Also, T xi 

Pi = 4a C 



4a C ib "+ 
+ C zW 



C ia C ia 



c ib C ib: 



and 



The constant A is the electron- 
phonon coupling related to the Jahn- Teller distortion of 
the MnC>6 octahedron^' 2 Regarding the phononic stiff- 
ness, and in units of t^ a — 1, the D a parameters are 
D\ = 1 and D2 = D3 = 0.5, as discussed in previous 
literature^ The rest of the notation is standard. In our 
effort here, the emphasis is on the case A=0 believed to 
be of sufficient relevance to deserve a special study since 
it already contains^ a competition between FM metallic 
and CE insulating states at (n)=0.5. Thus, this is an ex- 
cellent testing ground for the TPEM, particularly having 
in mind the next challenge involving a TPEM study in 
the presence of quenched disorder. However, briefly some 
results at nonzero A will also be shown. 
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FIG. 11: (color online). Convergence of the structure factors, 
varying the parameters of the TPEM algorithm: (a) Magne- 
tization \M\=y/S(0, 0) vs. T/t at J AF =0.0 using the DIAG 

r = 10~ 6 , and the values 



method and TPEM with e D 



=10" 



of M indicated, (b) Order parameter associated with the CE 
phase \O C e\ = ^S(tt,0) vs. T/t at Jap =0.2 using the DIAG 
method and TPEM with e pr = 10~ 5 , e tr = 10~ 6 , and the values 
of M indicated, (c) Magnetization \M\=^/S^0) vs. T/t at 
Jaf = 0.0 using the DIAG method and TPEM with M=20 
and etr = 10~ 6 , varying e pr —e as indicated, (d) Order param- 
eter of the CE phase |Oce|= \/S(tt, 0) vs. T/t at Jaf = 0.2 
using the DIAG method and TPEM with M=20, e tr = 10 -6 , 
varying e pl —e as indicated. All calculations were done on a 
12 x 12 lattice, using 1000 Monte Carlo steps for thermaliza- 
tion and 1000 for measurements. 



B. TPEM Performance 



Test of the TPEM in Small Systems 



As in the case of the one-band model, the analysis 
starts here by comparing DIAG and TPEM results on 
small systems. Figure ITTI contains the magnetization (in 
absolute value, and coming from the classical spins) vs. 
temperature. The results in (a) and (c) were obtained 
at Jaf=0.0, A=0, and (n)=0.5, a regime known to de- 
velop ferromagnetism at low temperatures^ Indeed, both 
methods show a nonzero value for the magnetization. 
The dependence with the TPEM parameters indicates 
that a M of approximately 30 or higher is sufficient to get 
accurate results. This is a conclusion that also appears 
in (b) where the case of a CE state is studied, which 
is stabilized with increasing Jaf- Regarding the other 
TPEM parameters, (d) shows that e pr =10~ 5 , as used for 
the one-band case, leads to accurate results. Overall, it 
seems that the same set of parameters deduced from the 
one-band model investigations can also be used for two 
bands, an interesting simplifying result. 
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FIG. 12: (color online). Lattice size dependence of the square 
root of the structure factors at the momenta characteristic 
of (a) a FM state, k=(0,0) and J AF =0.0, and (b) a CE 
phase, k=(7r,0) and Jaf=0.2. Results were obtained with 
the TPEM with M=20, e pr = 10" 5 , e tr =10" 6 . In addition, for 
(a) the DIAG method was also used on a 8 x 8 lattice as 
indicated. In the simulation, 1000 MC steps were used for 
thermalizations and 1000 steps for measurements. 



2. Dependence of Results on Lattice Sizes 

Figure ^| illustrates the dependence of results on lat- 
tice sizes. The systematic behavior is similar to that 
observed in the case of the one-band model at the same 
density (n)=0.5. In (a) results for the magnetization vs. 
temperature indicate the existence of a FM state at low 
temperatures, as well as small finite-size effects when the 
20x20 and 32x32 clusters are compared. Even less pro- 
nounced size effects are found in the CE regime, increas- 
ing Jaf as shown in (b). There is no indication that the 
TPEM deteriorates with increasing lattice size, providing 
hope that this method will be strong enough to handle 
the introduction of quenched disorder in future studies. 



C. Phase Diagram 

1. Results without Phonons 

To further test the TPEM, the phase diagram of the 
two-band model at A=0 and (n)=0.5 was obtained. Also 
at very low temperature, the energy was found as a func- 
tion of Jaf- The results are in Fig. Part (a) shows 
an excellent agreement among the several lattice sizes 
studied here. The abrupt change in the slope of the 
curve near Jaf=0.15 indicates a first-order transition, 
similar to that found in the one-band case and in previ- 
ous literature. In (b), the full phase diagram is obtained. 
There are clear qualitative similarities with the results 
presented before by Aliaga et al. using a 4x4 cluster £ In 
particular, the curve defining the CE phase at low tern- 
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FIG. 13: (color online), (a) Total Energy vs. Jaf at low 
temperature (T=0.01i) for different lattices as indicated. For 
the 12 x 12 lattice the DIAG method was used, while for 
the others the TPEM was employed with M=20, e pr = 10~ 5 , 
etr=10 -6 . In the simulation, 1000 MC steps were used for 
thermalizations and 1000 steps for measurements, (b) Phase 
diagram of Hamiltonian Eq. {y|j varying temperature and Jaf 
(A=0). The three magnetically different regions: FM, PM, 
and CE are indicated. The phase diagram was calculated for 
different lattices as shown. For 12 x 12 the DIAG method was 
used and for the others the TPEM with M=20, e pr =10 -5 , 
e t r = 10 -6 . The critical temperatures were obtained from the 
calculation of structure factors, as shown in Fig. 1121 
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FIG. 14: Phase diagram of Hamiltonian Eq. ® varying tem- 
perature and Jaf for A=0.5 and the harmonic parameters 
discussed in the text, i.e., with the inclusion of phonons. The 
critical temperatures were obtained from the calculation of 
structure factors on a 12 x 12 lattice with the DIAG method. 
Note the similarity of this phase diagram with the result ob- 
tained at A= 0.0, suggesting that to simulate the competition 
between FM metallic and CE insulating regimes the presence 
of a robust electron-phonon coupling is not necessary. 



2. Influence of Phonons 



perature has a positive slope rather than being vertical 
as in other cases. The fact that the TPEM gives results 
in excellent agreement with DIAG but on substantially 
larger systems is very encouraging and establishes this 
technique as a key method for a frontal attack to the 
CMR problem using realistic models and quenched dis- 
order. 

For completeness, CPU times for the case of the two- 
band model are provided in Table ITT1 The CPU time per 
site does not change dramatically with N, close to the 
expected theoretical estimation for the TPEM. Clearly, 
lattices well beyond 32 x 32 can be handled with this tech- 
nique. 



TABLE II: CPU Times of the TPEM in seconds per 5 Monte 
Carlo steps for Hamiltonian Eq. ||]|) with Jaf=0 and inverse 
temperature j3 = 50, and the lattices shown. The third col- 
umn is the ratio of CPU time per lattice site. The computer 
used was an AMD Opteron(tm) 244, 1.8GHz with 1MB cache. 
The TPEM parameters were M=20, e pr =10~ 5 , and e tr = 10 -6 . 



LxL 


CPU Time(s) 


CPU Time/N 


12 x 12 


124 


0.86 


20 x 20 


642 


1.61 


24 x 24 


992 


1.72 


32 x 32 


2451 


2.39 



As explained in the introduction, the general scenario 
proposed for manganites does not depend on particular 
details of the competing phases, but in the competition 
itself. 2 Thus, to the extent that the FM metallic and CE 
insulating phases are found in competition, the value of 
electron-phonon coupling A is not of crucial relevance. 
We believe that A is likely small in practice, since recent 
experiments are not finding evidence of a robust charge 
checkerboard (see, for example, Ref.l25h and, in addition, 
theoretical studies have shown that a large A renders 
the FM state also insulating^ Nevertheless, to confirm 
that the results are not severely affected by switching on 
A, in Fig. ^] the phase diagram for A=0.5 is presented 
on a lattice substantially larger than used in previous 
investigations. 9 Comparing Figs. 1131 and !14l clearly both 
cases lead to very similar phase diagrams. Since remov- 
ing the phononic degree of freedom speeds up the simu- 
lations, these results suggest that the future effort in this 
context could focus in the A=0 case and still expect to 
find realistic conclusions. 



D. Density of States 

As in the case of the one-band model, we also tested 
whether the TPEM technique can reproduce the DOS of 
the two-band model in the regime where the system is 
insulating (CE phase). The results are in Fig. Part 
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FIG. 15: (color online) (a) DOS of a perfect CE phase (sin- 
gle spin configuration) obtained on a 12x12 lattice calculated 
with the DIAG method (solid black line) and with the TPEM 
(dashed red line) using M=100. The chemical potential lies 
in the left gap (arrows) indicating that the system is an in- 
sulator, (b) DOS of the system with Jaf=0.2 (CE-phase 
ground state) on a 20x20 lattice calculated with the TPEM 
using M=100, as described in the text. The location of the 
chemical potential is indicated by the vertical dashed line. 



(a) shows a comparison between DIAG and TPEM on a 
12x12 cluster. The agreement is excellent for the case 
of M=100 (shown), and fairly acceptable for smaller val- 
ues of M. For larger lattices that can only be studied 
with TPEM (part (b)), the results are also in good agree- 
ment with expectations. Then, no problems have been 
detected in calculating the DOS using the TPEM tech- 
nique in the regime where the model is in an insulating 
state. 



E. Conductances: Comparison TPEM vs. DIAG, 
and Results for Increasing Lattice Sizes 

As in the case of the one-band model, the final test for 
the two-band case is the calculation of the conductance. 
Results are shown in Fig. 1161 Even using clusters much 
larger than can be handled with DIAG, the behavior of 
the conductance is properly captured by TPEM. There 
are no hidden subtleties involved in this estimation of 
transport properties, opening the way toward many fu- 
ture applications of this technique. 



F. Influence of Magnetic Fields in the Clean Limit 

Similarly as for the case of the one-band model, we 
have also studied the influence of a magnetic field in the 
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FIG. 16: (color online) Upper Panel: Conductance vs. tem- 
perature for Jap = 0.0 (ferromagnet) and Jaf = 0.2 (CE 
phase), different lattices sizes and algorithms, as indicated. 
Lower Panel: Logarithm of the resistivity vs. temperature 
for the same parameters as before. 



case of the two-band model. The value of Jaf was chosen 
to be 0.175, namely on the CE side but close to the first- 
order metal-insulator transition. As shown in Fig. El 
the application of a relatively small field - in the natu- 
ral units of the problem - leads to a drastic change in 
the resistivity at low temperatures. In this regime, the 
insulator is transformed into a metal (negative magne- 
toresistance) . As remarked before, this is in agreement 
with previous studies carried out by Aliaga et aim, show- 
ing that the "low temperature" large magnetoresistance 
materials^ can be understood by double-exchange mod- 
els in the clean limit. However, as in the case of one-band 
models studied before, the finite-temperature CMR ef- 
fect in Mn-oxides cannot be understood with clean limit 
models, as shown in Fig. El since the zero magnetic-field 
resistivity does not have the large peak characteristic of 
manganites. Future work with TPEM will address this 
issue including quenched disorder. 



V. CONCLUSIONS 

In this paper, it has been shown that the computa- 
tionally intensive exact-diagonalization algorithm for the 
study of CMR-manganite models can in practice be re- 
placed by the novel Truncated Polynomial Expansion 
Method (TPEM) without any substantial loss in accu- 
racy. The DIAG, although being exact, does not permit 
simulations of large clusters owing to the fact that the 
computational cost grows like the 4 th power of the sys- 
tem size, N. On the other hand, the newly developed 
TPEM algorithm reduces the computational complexity 
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FIG. 17: (color online) Effect of a small magnetic field on 
the resistivity of the CE phase, for couplings in the vicinity 
of the first order transition metal-insulator. Shown is the 
resistivity vs. T at Jaf =0.175, A=0, on a 20x20 lattice for 
B=0.1t and without field (B=0) for comparison. Note the 
large change in the resistivity at low temperature, compatible 
with a colossal magnetoresistance. The TPEM was used with 
M=30, e pr =10~ 5 and e tr = 10" 6 . 1,000 Monte Carlo steps 
were done for thermalization and an additional 100 more for 
measurements. 



to 0(N), thereby allowing for simulations of larger sys- 
tems. 

For both the one- and two-band double exchange 
model with interacting Mn spins, we have compared sys- 
tematically the results calculated with both the diago- 
nalization method and the TPEM algorithm. As the 
spin-spin coupling Jaf is varied, the one-band model re- 
veals a low-temperature first-order phase transition be- 
tween conducting FM and insulating Flux states in the 
vicinity of Jaf=0.045. For two bands, a similar first- 
order transition separates FM metallic and CE insulat- 
ing phases. Our calculations presented here included a 
systematic study of the performance of the TPEM al- 
gorithm varying its parameters M, e pr and et r , already 
defined in the introductory sections. It was shown that 
the results of the TPEM algorithm converge to those of 
the DIAG algorithm with increasing M, and for M > 30 
TPEM results are sufficiently accurate to obtain phase 
diagrams with small error bars. This number appears to 
be fairly stable under changes in the model, couplings, 
and for different phases, and lattice sizes. Also nothing 
indicates that working at density different from 0.5, the 
focus of the current effort, will spoil the TPEM perfor- 
mance. [However, it is advisable to be particularly cau- 
tious near critical temperatures, where in some cases we 
found the need to increase M to 40] . Similar systematic 
results were presented for the e's. Overall, the general 
process of fixing TPEM parameters on small systems by 
comparing with DIAG and then using the same parame- 



ters on larger lattices appears reliable, and this method 
will be applied to other systems in the near future. 

Parallelization of the TPEM algorithm makes it pos- 
sible to study large clusters. Taking advantage of this 
possibility, in the absence of quenched disorder, we have 
made calculations on lattices of up to 40 x 40 sites even for 
the case of a finite Hund coupling in the one-band model. 
But previous calculations in the limit of Jh=oo used up 
to 8000 sitesj2& thus even accounting for the factors of 2 
involved in comparing finite and infinite Jh there is still 
room for further improvement. We have also checked that 
calculations of conductances, crucial to predict transport 
properties, also can be carried out smoothly with the 
TPEM, and in addition the size effects are in general 
small. While obviously the increase in lattice size al- 
lows for more accurate determinations of critical temper- 
atures, even more importantly these large lattices will be 
crucial for the next big step in large-scale manganite sim- 
ulations which is the introduction of quenched disorder. 
This disorder is expected to lead to a percolative-like pic- 
ture that causes the CMR phenomenon. Any percolative 
effect requires large systems, and having access to clus- 
ters substantially larger than those studied with DIAG is 
a key necessary condition to unveil the conceptual reason 
behind the CMR phenomenon. 

Interesting physical results were also here reported. 
This includes the one-band phase diagram with a metal- 
insulator transition. But the main result in this context is 
the presence of an enormous magnetoresistance effect at 
low temperatures, even in the clean limit studied here. 
This effect was already anticipated by Aliaga et al. in 
their pioneering work on this subject on small systems^ 
The survival of the effect on the large clusters reachable 
by the TPEM, as described in this manuscript, shows 
that some forms of CMR found experimentally can al- 
ready be accurately reproduced using realistic models, pro- 
viding further support that a theoretical solution of the 
CMR puzzle is within reach. But for the most common 
form of CMR in Mn-oxides at temperatures close to the 
Curie temperature, studies with quenched disorder will 
likely be needed. 

Summarizing, we here reported a successful implemen- 
tation of the TPEM for the study of double-exchange-like 
models for manganites. The technique has a CPU time 
that grows linearly with the number of sites, and in ad- 
dition it is parallelizablc. Thus, the main result is that a 
novel technique has been identified and tested that can 
led to a frontal attack of the most interesting problem in 
manganites: the analysis of large magnetoresistance ef- 
fects in the presence of quenched disorder, when phases 
compete via a first-order transition in the clean limit. 
This "holy grail" of simulations will be the focus of our 
effort in the near future. It will demand at least an order 
of magnitude more effort than in the present manuscript, 
but this can be alleviated by increasing the number of 
nodes available for the simulations and we are already 
working on this aspect. The large-scale computational 
facilities at Oak Ridge National Laboratory will play a 
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key role in reaching this ambitious goal. 
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